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ABSTRACT 

We discuss the link between dark matter halos hosting the first PopIII stars 
and the rare, massive, halos that are generally considered to host bright quasars 
at high redshift {z ^ 6). The main question that we intend to answer is whether 
the super-massive black holes powering these QSOs grew out from the seeds 
planted by the first intermediate massive black holes created in the universe. 
This question involves a dynamical range of 10^^ in mass and we address it 
by combining N-body simulations of structure formation to identify the most 
massive halos at z ~ 6 with a Monte Carlo method based on linear theory to 
obtain the location and formation times of the first light halos within the whole 
simulation box. We show that the descendants of the first ~ IO^Mq virialized 
halos do not, on average, end up in the most massive halos at z 6, but rather 
live in a large variety of environments. The oldest PopIII progenitors of the most 
massive halos at 2; ~ 6, form instead from density peaks that are on average one 
and a half standard deviations more common than the first PopIII star formed 
in the volume occupied by one bright high-z QSO. The intermediate mass black 
hole seeds planted by the very first PopIII stars at 2; > 40 can easily grow to 
masses msH > lO^'^M© hj z = 6 assuming Eddington accretion with radiative 
efficiency e < 0.1. Quenching of the black hole accretion is therefore crucial to 
avoid an overabundance of supermassive black holes at lower redshift. This can 
be obtained if the mass accretion is limited to a fraction rj ^ 6 ■ 10^^ of the total 
baryon mass of the halo hosting the black hole. The resulting high end slope of 
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the black hole mass function at z = 6 is a ~ —3.7, a value within the la error 
bar for the bright end slope of the observed quasar luminosity function at z = 6. 

Subject headings: cosmology: theory - galaxies: high-redshift - early universe - 
methods: N-body simulations 



Introduction 



Bright quasars at z Q are very luminous and rare objects that can can be detected out 
to huge cosmolo gical distances in very large area surveys like the Sloan Di gital Sky Surve y 
(IPan et all 12004 1 Their estimated space density is ^ 2.2 • lQ-'^{Mpc/h)-^ flPan et allbooj ). 
that is about one object per about 200 deg^ of sky , assuming a depth of = 1 centered at 
z = 6.1 under the third year WMAP cosmology (ISpergel et al.l 120061 ) . Their luminosity i s 
thought to be due to accretion onto a super-massive black hole (e.g., see lHopkins et al.ll2005l ). 
A common expectation is that the luminous high-z quasars sit at the center of the biggest 
proto-clusters at that time. Some observational evidence of over-densiti es of galaxies in two 



deep HST-ACS f ields containing a bright z=6 quasar has been claimed (jStiavelli et al.l 12005 



Zhen et al.l 120061 ). but it is unclear whether this is true in general. In fact an ACS image 
only probes a long and narrow field of view of about 6 x 6 x 320(Mpc//i)^ in the redshift 
range [5.6 : 6.6], so a significant number of detections may come from galaxies unrelated to 
the environment of the host halo of the bright quasar. 

Numerical simulations to address the formation of bright quasars are extremely chal- 
lenging given their low number density. A huge simulation cube with edge of > 700 Mpc/h 
is required just to expect, on average, one such object in the simulat ion box. A major com- 
putational investment, like the Millennium run (ISpringel et al.ll2005l ). is required to resolve 
at high redshift {z > 20) virialized halos on this volume and to follow their merging history 
down to ^ ~ 6. Even assuming that the simulation volume is big enough that there is the 
expectation to find halos hosting bright quasars, how can these halos be identified? In prin- 
ciple two, non mutually exclusive, alternatives appear plausible: either the super-massive 
black holes are hosted in the most massive halos with the corresponding number density of 
SDSS quasars or these black holes have grown from the first PopIII Intermediate Mass Black 
Hole seeds, therefore representing the descendant of the rarest density peaks that hosted 
first stars. 



The first scenario implies that the rri BH—o' relation 



is already in place at high redshift (IVolonteri et al.l 



Perrarese fc Merritt 



2003 



Hopkins et al. 



200G 



2005 



Gebhardt et al. 



Di Matteo et al, 



2005l ). In that case multigrid simulations can be carried out to follow in detail the growth 
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of the supermassive black hole (e.g., see iLi et al.ll2006l ). In the second scenario the quasars 
progenitors would be traced back to the first PopIII stars created in the universe within 
«i IQ^Mq mass halos virialized at z 50 (jSromm &: Larson 2004 : Abel et al. 2002). These 
PopIII stars are very massive M > IOOMq, so after a short life of a few million years explode 
and may leave intermediate mass black holes, plausible seeds for the super-massive black 
holes observed at lower redshift. Of course the two scenarios can be consistent with each 
other if the first p erturbations to collap se are also the most massive at z = Q. This seems to 
be implied, e.g. in lSpringel et al.l (120051 ) . where the bright quasar candidate in the simulation 
is traced back to one of the 18 collapsed halos at z = 16.7. 

In this paper we explore the link between the first PopIII halos collapsed in a simulation 
box and the most massive halos at lower redshifts to gain insight on the scenarios of bright 
quasar formation. This is a numerically challenging problem as the dynamical range of 
masses involved is very large: a simulation volume of 5 • 10^{Mpc/h)^ has a mass of about 
3.3 • IO^^Mq/Zi, that is more than 10^^ times the mass of a PopIII dark matter halo. We 
have adopted an original approach to the problem, broadly inspired by the tree method by 
Cole et al.l (1l994j ). We first simulate at relatively low resolution the evolution of a simulation 
volume down to z = 0. Then, starting from the density fiuctuations field in the initial 
conditions of the numerical simulation, we compute analytically the redshift distribution of 
the oldest PopIII halo collapsed within each single grid cell. The information is then used 
as input for a Monte Carlo code to sample for each particle of the simulation the collapse 
redshift of the first PopIII progenitor dark matter halo. The formation time of the oldest 
PopIII remnant within the most massive halos identified at z ^ 6 is finally compared with 
that of the oldest PopIII star sampled over the whole simulation volume and the implications 
for the growth of supermassive black holes are discussed. 

Our approach is tuned to investigate the formation and the subsequent remnant distri- 
bution of the first, rare density peaks that hosted PopIII stars at ^ > 30. With this respect 
our study has a similar goal to iReed et al.l (120051 ) . with the important difference that we 
search for the first PopIII star in the complete simulation box and not by means of pro- 
gressive refinements around substructures that probe only a small fraction of the total box 
volume. As our method is tuned at finding very rare fiuctuations, it is not easily applied to 
the significantly more common 3a peaks with mass ~ IO^Mq that collapse at z ~ 20 and 
that might constitute the majority of PopIII sta rs, if these are terminated by chemical feed- 
back at z < 20 (e.g. see. iGreif fc Bromml 120061 ) and not by photo-dissociation of molecular 
hydrogen at 2 > 25 (IHaiman et al.ll2000l ). 



This paper is organized as follows. In Sec. [2] we present the details of the numerical 
simulations carried out. In Sec. [3] we analyze the numerical results focusing on the merging 
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history of the first PopIII halos formed in the simulation box. In SecJUwe review when the 
first stars epoch end, while in Sec. [5] we discuss the imphcations of the PopIII distribution 
that we find for the build-up of supermassive black hole population at 2; < 6. We conclude 
in Sec. Ei 



2. Numerical Methods 

2.1. N-body simulations 

The numerical simulations presented in t his paper have been carried out using the 
public version of the PM-Tree code Gadget-2 (ISpringell 120051). Our standard choice is to 
adopt a cosmology based on the third year WMAP data (ISpergel et al.ll2006l ): fl^ = 0.74, 
flrn = 0.26, -f^o = 70 km/s/Mpc, where Q^n is the total matter density in units of the 
critical density (pc = ^Hq/{87tG)) with Hq being the Hubble const ant (paramete rized as 
Ho = lOOh km/ s/Mpc) and G the Newton's gravitational constant (jPeebled Il993l ). is 
the dark energy density. As for erg, the root mean squared mass fluctuation in a sphere of 
radius 8Mpc/h extrapolated at 2; = using linear theory, we consider both cxg = 0.9 and 
(Tg = 0.75, focusing in particular on the higher value that p rovides a better match to the 
observed clustering properties of galaxies (lEvrard et al.l 120071 ). 



The initial con ditions have been generated using a code based on the Grafic algorithm 
( iBertschingerl I2OOII ). An initial uniform lattice is perturbed using a discrete realization of 
a Gaussian random field sampled in real space and then convolved in Fouri er space with a 
ACDM transfer function computed using the fit by lEisenstein fc Hul (119991 ) and assuming 
a scale invariant long- wave spectral index (n = 1). The initial density field is saved for 
later reprocessing through the first light Monte Carlo code (see Sec. 12.21) . The particles 
velocities and displacements are then evolved to the desired starting redshift {zstart = 65.67, 
i.e. ttstart = 0.015) using the Zel'dovich approximation and the evolution is followed using 
Gadget-2 (ISpringell l2005l ). Dark matter halos are identified in the simulations snapshots 
using the HOP halo finder (lEisenstein &: Hutlll998l ). 



To find the optimal trade-off between mass resolution and box size, both critical param- 
eters to establish a connection between PopIII halos and the most massive halos identified 
at 2; = 6, we resort to simulations (see Tab. [1]) with three different box sizes, all simulated 
with = 512^ particles: 



(i) A "large" box size of edge 720 Mpc/h that is large enough to contain on average about 
one bright high-z quasar. The mass resolution is 3.7 • IO^^Mq/Zi (corresponding to a 
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halo of 20 particles). 

(ii) A "medium" box size of edge 512 Mpc/h that represents a compromise between a 
slightly higher mass resolution than (z) and a still reasonably large simulation volume. 

(iii) A "small" box size of edge 60 Mpc/h. While this box size is too small to host a 
bright z K, Q quasar, i ts volume is still larger than that of deep surveys like the UDF 



( iBeckwith et al.ll2006l ). that spans a volume about 20 times smaller than this box in 
the redshift interval z e [5.6 : 6.6] (the typical redshift uncertainty for i-dropouts). 
Halos down to about 3 • IO^Mq/Zi can be identified in this box. The analysis of the 
results from this simulation will show the fundamental role played by the large volume 
employed for simulations (i) and (ii). 



2.2. Monte Carlo code for first light sources 



Given the initial density fluctuations field on the simulation grid, where a cell has a 
mass of order 10^'^ — 10^^Mq//i, our goal is to estimate the redshift of the first virialized 
perturbation wi thin each cell at the ni ass scale of early PopIII dark matter halos (i.e. < 
lO^M©, see e.g. iBromm fc Larsorul2004j ). For this we resort to an analytical treatment based 
on a linear approximation for structure formation. 

The initial conditions for a N-body simulation in a box of size L with particles and a 
single particle mass mp define a Gaussian random field 5p for the N cells (associated to the 
location of the N particles) of the simulation grid. This density field is usually generated by 
convolving white noise with th e transfer function associated to the power spectrum of the 
density perturbations (e.g., see iBertschingerl 120011 ) and is used to obta in the initial v elocity 
and positions displacements for the particles (e.g. see Eq. 5.115 in iPeebled Il993l ). The 



density fluctuation in each cell has a contribution from different uncorrelated frequencies in 
the power spectrum. When the initial conditions for an N-body simulation are generated, 
the power spectrum has an upper cutoff around the Nyquist frequency for the grid used (i.e. 
around the frequency associated to the average inter-particle distance) and a lower cutoff at 
the frequency associated to the box size (if periodic boundary conditions are enforced). A 
higher resolution version of the initial density field can be obtained by simply increasing the 
grid size and adding the density perturbations associated to the power spectrum between 
the old and the new cutoff frequencies. 

In linear approximation one can use the field 5p to obtain the redshift of virialization of 
a structure of mass Mh > rrip at a given position x in the grid. To do this one averages the 
field 6p using a spherical window centered at x with a radius such that the window encloses 
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a mass Mh and computes assuming linear growth the redshift at which the average density 
within the window reaches 6p = 1.69 (in units where the average density of the box is 1). 
In fact, for a spherical collapse model, when 6p = 1.69 in linear theory, then the halo has 
reached virial equilibrium under the full non-linear dynamics. This concept is at the base of 
the various pr oposed methods for computing analytically the mass function of d ark matter 



halos (e.g., see lPress fc Schechterlll974l : iBond et al.lll991uSheth fc Tormenlll999l ). 



We apply this idea to estimate the formation rate and the location in the simulation 
volume of dark matter halos at a mass scale below the single particle mass used in the 
simulation. A straightforward implementation consists in generating first the density field 
associated to the N-body simulation, and then to refine at higher resolution the field by 
means of a constra ined realization of the initial conditions used in the N-body run (e.g. see 



Bertschingerll200ll ) . This provides exact and complete information on the whole density field, 
but the price to pay is the execution of very large Fast Fourier Transforms on the refined 
grid. If the goal is to compute density fluctuations down to a mass of ~ 10^ Mq/H over a box 
of edge 720Mpc/h, a grid of 29184^ is needed, which would require about 100 TB of RAM, 
that is well beyond the current memory capabilities of the largest supercomputers. 

A shortcut is however available, if one trades information for numerical complexity. 
Given a realized numerical simulation, we are in fact not interested in getting a detailed 
picture of the dynamics at sub-grid resolution, but only in identifying for each grid point the 
redshift of virialization of its first progenitor at a given sub-grid mass scale. For example, 
given a simulation with single particle mass of 10^°Mq//i our aim is to quantify the redshift 
of virialization of the first dark matter halo of mass IO^Mq/H within the volume associated 
to the IO^^Mq/H particle. In that case, if we were to have the full sub-grid information we 
would search for the maximum realized value of the density within the 10^ sub-grid cells of 
mass 10^ Mq/H that constitute our IO^^Mq/Zi single particle cell. As the density fluctuation 
field is a Gaussian random field, the density in sub-grid cells will be a Gaussian centered at 
the density of the parent cell and with variance given by integration of the power spectrum 
of density fluctuations truncated between the Niquist frequency of the parent cell and that 
of the sub-cells. 

Therefore, for a single cell of mass rrip, the redshift of collapse of a sub-grid progenitor 
at a mass scale rrifs can be obtained simply by sampling from the probability distribution of 
the maximum of the sub-grid fluctuations of the k = mp/rrifs sub-cells of mass rufs that are 
within a cell of mass rup. The probability distribution for the maximum of these fluctuations 
is available in analytic form when the field is Gaussian, as in the case considered here. In 
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fact, given a probability distribution p{x), with partition function P{x): 

/X 
p{a)da, (1) 
-oo 

the probabihty distribution q{m, k) for the maximum m of k random numbers extracted 
from p{x) (m = Max{xi, ...,Xk)) is the derivative of the partition function for m, that in 
turns is simply the k-th power of the partition function for ). Therefore we have: 

q{m,k) = k ■ p{m) ■ P{mf'^. (2) 

Eq. ([2]) has a simple interpretation: the probability that the maximum of k random numbers 
lies in the interval [m, m + Sm] is given by the probability of sampling one of the k numbers 
exactly in that interval and all the other numbers below m. 

With the aid of Eq. ([2]) we can sample the distribution of the maximum of the additional 
sub-grid density fluctuations that need to be considered in order to probe the mass scale 
of PopIII halos. The variance cx/s to be used in p{x) may be computed from the power 
spectrum of the density fluctuations by considering an upper cut-off at the wavelength of 
one cell size in the initial conditions. Or, equivalently, if the complete power spectrum of 
density fluctuations has variance (T{Mgrid_ceii) at the mass scale of one grid cell and variance 
cr{Mpopin_haio) at the mass scale of a halo hosting a first star, we set ajs such that: 

= Or{MpopIII_haloY — CT{Mgrid_ceuY ■ (3) 

Therefore our recipe for estimating the age of the earliest progenitor formed in each cell 
is the following: 

(i) Starting from the initial density fluctuation field on the grid used to initialize the 
N-body run compute the mass refinement factor k to go from the mass of a sin- 
gle particle (i.e. the mass within one grid cell) to that of a PopIII star halo {k = 

Mgrid_cell/ M PopIII Jialo ) ■ 

(ii) Given the power spectrum of the density fluctuations, Mgrid_ceii and Mpopiii_haio com- 
pute ajs- 

(iii) Extract one random number r from q{m, k) (see Eq. [2]) where p{m) is a Gaussian 
distribution with zero mean and variance afs- 

(iv) Sum r to the value of the density field in the cell to obtain (5p/s/ p)max in the cell. From 
the value of {6pfs/ p)max it is then straightforward to compute the non-linear redshift 
for that perturbation, i.e. the redshift Zni when the linear density contrast reaches a 
value {6pfs/p) max (Znl) = 1.69. 
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The particles of the simulation now carry the additional information of the redshift at 
which their first PopIII star dark matter halo progenitor has collapsed in linear theory (a 
proxy for the redshift of actual virialization). Once halos have been identified in simulation 
snapshots, the redshift of the earliest PopIII progenitor within the halo is easily obtained. 
It is similarly easy to identify in a snapshot what is the environment in which the particles 
with the oldest progenitors live. This procedure is robust with respect to variations of 
the simulation resolution, as long as the focus is on rare density peaks, with an average 
occupation number per simulation cell (i.e. particle) much smaller than unity. Numerical 
tests are presented in Appendix Rl 

This method has two main advantages: 



It allows to use relatively inexpensive "low resolution" simulations to identify the 
largest objects at low redshift {z < 6). In fact if we are interested in identifying 
the most massive halos at 2; ~ 6 as host halos for quasa r candidates a mass resolution 



of ~ 10 Mq/Zi is sufficient (e.g. in lSpringel et al.ll2005l the mass of the largest halo at 



z = 6.2 is 3.9 • W'^Mq/H for a simulation volume of {500Mpc/hf). 

(ii) For a given numerical simulation, several Monte Carlo realization can be generated to 
gather robust statistical constraints on the properties of dark matter halos hosting first 
light sources as well as the spatial distributions of the first halo remnants in halos at 
lower redshift. 



However our method has the drawback that it cannot be easily extended to the inves- 
tigation of the detailed merging history at the sub-grid level, as only the virialization time 
of the earliest progenitor of each particle at a given mass scale is provided. In addition, 
the identification of the first virialized PopIII halos is expressed in terms of the halos with 
the highest Zni- We are therefore neglecting the non-linear evolution and the environmental 
dependences on the dynamics of the dark matter collapse, such as tidal forces, therefore 
missing the precise redshift at which a PopIII halo virializes. These are limitations that we 
need to accept as the non linear evolution could be followed over the whole box only at the 
price of running a simulation prohibitively intensive in cpu and memory res ources, with at 
least 10^ time more particles than in the Millenium Run (jSpringel et al.ll2005l ). This appears 
unfeasible for the time being, ev en considering next generation dedicated supercomputers, 
like the GrapeDR jMakinolboOsh . 
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3. The fate of the first PopIII halos 

3.1. Analytical considerations 

The general picture for the connection between first halos and the most massive halos 
at z ~ 6 can be obtained using analytical considerations, that will be later confirmed in 
Sec. 13.21 by the results of our numerical investigation. 

Following the choice for our large box simulation, we consider a volume of {720Mpc/h)^ 
of mass Mbox, large enough to host a bright z ^ 6 quasar. We estimate from the Press- 
Schechter formalism (see also the masses of the z = 6 halos in our "large" box simulation 
in Sec. 13. 2p that the most massive halo at z ~ 6 has a mass (that we call Mqh) of about 



10 — 10 Mq/Zi (see also ISpringel et al.l l2005l ) . Since the most massive halo is the first at 



its mass scale to be formed, through the use of Eq. [2] we can obtain the distribution of its 
initial density fiuct uation (see Fig . [T])- If we assume Mg^ = = 4.3 • IO^^Mq/Zi 



(in agreement with ISpringel et al.l 120051 ). this halo is expected to have originated from a 
density fiuctuation in the range [5 : 5.7]a{Mqh) at 90% of confidence level. We now consider 
the volume initially occupied by the mass Mgh and we compute from the primordial power 
spectrum the variance afg of density perturbations at mass scale of a Poplll halo [Mfg = 
IO^Mq/H = 1/160^ Mqh) considering only contributions from wavelengths at a scale below 
the volume enclosed by Mgh (see Eq. [3]). We obtain ct/s = 4:.85a{Mqh). From Eq. [2] follows 
that the maximum of 160^ Gaussian random numbers with variance ajs is distributed in the 
range [23.4 : 27.0]a{Mqh) at 90% of confidence level. Combining the two 90% confidence 
level intervals, this means that the first PoplII progenitor of a bright quasar originated from 
a perturbation in the range [28.4 : 32.7]a{Mqh). If we consider instead a random sub-cell 
among the 160^, the probability that the maximum sub-grid perturbation is smaller than 
32.7a{Mqh) is only 0.99995, so several hundreds of the 180^ cells among the whole simulation 
volume are expected to have a PopIII progenitor formed before that of the most massive 
z = 6 halo. In fact from integration of Eq. [21 the sigma peak associated to the first star 
in the box is expected to be greater than 35.5a{Mqh) at 99.99 % of confidence level (and 
in the interval [36.2 : 38.8]a{Mqh) at 90% of confidence level). Therefore the rarity of the 
earliest PopIII progenitor of the most massive halo at z = 6 is about 1.50"/^ less than that 
of the first PopIII star formed in the simulation volume. In terms of formation redshift, the 
first PopIII star dark matter halo in the simulation volume virializes in the redshift interval 
2 G [49 : 53], while the earliest PopIII progenitor of the QSO halo is formed at z G [38 : 44] 
(both intervals at 90% of confidence level and computed for as = 0.9). 

The picture changes quite significantly if we consider a smaller box size. E.g. in our 
SI simulation (see Tab. [1]) with a volume of (GOMpc/Zi)^ a perturbation on a mass scale 
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Mqh = 7 • IO^^Mq/H ~ Mbox/'27^ is expected to be the most massive at z ^ 6. Such a 
halo derives at 90% of confidence level from a fiuctuation [3.6 : 4.6](j{Mqh). If we further 
assume Mfg = S.S-IO^Mq/Zi = l/94^Mq/j, we have a sub-grid variance ajs = 3.66 so that the 
maximum of the 94^ random first light perturbation in a cell of mass Mgh is distributed in 
the range [16.4 : 19.4:]a{Mqh) at 90% of confidence level. By combining the two intervals as 
above, we expect that the first PopIII progenitor of the most massive z = 6 halo derives from 
a [20.0 : 2A.Q]a{Mgh) peak. The first PopIII star derives instead from a [23.8 : 2Q.l]a{Mgh) 
peak (always 90% of confidence level). At variance with the larger box, here the correlation 
between the most massive halo at 2; = 6 and the first PopIII star in the simulation is expected 
to be stronger and the most massive halo is likely to have as progenitor one of the first 10-100 
Pop III stars. 

From these simple analytical estimates it is clear that the most massive and rarest 
structures collapsed around z ~ 6 do not descend from the rarest sigma peaks at the first 
light mass scale in the simulation volume, when the simulation box represents a significant 
fraction of the Hubble volume. Conversely the black holes remnants of the first PopIII 
stars in the universe do not provide the seeds for super-massive black holes within the most 
massive halos at 2; < 6. The descendants of first PopIII stars are instead expected to be 
found at the center of a variety of halos, as we quantify in the next Section by means of 
N-body simulations. 



3.2. Simulations Results 

In cons tructing the z = Q halo catalogs we adopt the following parameters for the HOP 



halo finder (lEisenstein fc Hutlll998l ). The local density around each particle is constructed 
using a 16 particles smoothing kernel. For the regrouping algorithm we use: Speak = 240, 
^saddle = 180, 6outer = 100 and a minimum group size of 20 particles. In the large simulation 
box (run LI in Tab [1]) we identify 47 halos with 20 particles or more and the most massive 
halo (37 particles) has a mass of 6.9 ■ IO^'^Mq/H. In the medium simulation box (run Ml in 
Tab [1]) the higher mass resolution allows us to identify 694 halos with at least 20 particles 
and the most massive halo has 92 particles for a total mass of 6.1 • IO^^Mq/H, consistent 
with the results from the larger box. Finally in the small box simulations (runs SI and 
5*2 in Tab [T]) there are 14972 halos with at least 100 particles in 5*1 (cxg = 0.9) and 7531 
halos with at least 100 particles in 5*2 (erg = 0.75). The most massive halo has a mass of 
2.4 ■ 10^'^ Mq/H in 5*1 and of 7.1 • IO^^Mq/Zi in 5*2. The z = Q halo mass dist ribution for these 



two si mulations is well described (with displacements within ^ 15%) by a lSheth &: Tormen 



( 119991 ) mass function. 
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The link between the halos identified in the snapshots and the first hght sources is 
estabhshed using the Monte Carlo method described in Sec. I2.2[ For the large box we 
consider a refinement factor = 57^ to move from the single particle mass of the simulation 
to a typical PopIII halo mass, so that Mfg = 1.0 ■ 10^ M^/h. For the first 10 most massive 
halos at 2; = 6 we show in Fig. [2] the distribution of the redshift at which the oldest progenitor 
crosses the virialization density contrast threshold in linear theory {Sp/p = 1.69) and the 
distribution of the ranking of the collapse time computed over all the PopIII progenitors of 
the simulation particles. The collapse rank of the first PopIII progenitor of the most massive 
z = 6 halo is in the interval [474 : 45075] at 90% of confidence level, with median 8535. The 
corresponding virialization redshifts are in the interval [39.3 : 44.1] with median 41.1. For 
comparison the first PopIII halo in the box virializes in the redshift interval [49.0 : 52.7] with 
median 50.3; the 100th first light in the box collapses in the redshift range [45.5 : 45.8] with 
median 45.7. These results from the combined N-body simulation and Monte Carlo code 
are in excellent quantitative agreement with the analytical estimates of Sec. 13.11 and confirm 
that in a large simulation box the most massive halos at 2; = 6 do not derive from the rarest 
sigma peaks at the first light mass scale. This result is robust with respect to the adopted 
typical mass for PopIII halos. In Fig. Owe show the results obtained by changing the mass of 
the halos hosting the first stars considering larger halos {Mfs = 3.4 ■ IO^Mq/h with k = 38^) 
and smaller halos {Mfs = 3.0 ■ IO^Mq/h with k = 85^). The formation redshift varies as the 
first halos are formed earlier when they are less massive, but the relative ranking between 
the first PopIII halo in the box and the first PopIII progenitor of the most massive structures 
at 2; = 6 remains similar. In passing we note that our distribution of the formation redshift 
for the first 3 ■ IO^Mq progen itor of the most in assive z = 6 halo (formed at z ~ 46) is in 



agreement with the results by lReed et al.l (120051 ) . obtained by means of N-body simulations 



with adaptive refinements. However this halo is not the first one formed in the simulation 
box as we find that the first structure on this mass scale is formed at 2; > 55 (see Fig. [3]). 

The results are similar for the medium box, which has a volume that is only three 
times smaller than the large one (see Fig. Hj). The refinement factor used here is = 40 
that gives Mfs = 1.0 • 10^ Mq/H. The first PopIII halo in the box virializes in the redshift 
range [48.3 : 52.4] with median 49.7, while the oldest PopIII progenitor of the most massive 
halo virializes in the redshift range [38.5 : 43.2] (median 40.4) and has a collapse ranking in 
[574 : 37499] at 90% of confidence level, with median 7261. 

The picture changes significantly (see Fig. [5]) for the small box that has a volume more 
than 10^ times smaller than the large box. Here we use a refinement factor k = 5, that 
leads to Mfs = 8.6 • 10^Mq//i. The collapse rank of the first light progenitor of the most 
massive z = Q halo is in the range [1 : 103] at the 90% confidence level with median 13. The 
correlation between the first PopIII star and the most massive structure at z = 6 is therefore 
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strong due to the small volume of the box. This means that locally the oldest remnants of 
first stars are expected to be within the largest collapsed structures. 

From the medium box size numerical simulation we have also characterized the fraction 
of first PopIII remnants that end up in identified halos at 2; = 0. If we consider one of 
the first 100 first light halos collapsed in the box, there is an average probability of 0.72 of 
finding its remnant in a halo identified at 2; = with more than 100 particles (that is of 
mass above 6.7 • IO^^Mq/Zi). The median distribution for the mass of a halo hosting one of 
the remnants of these first light sources is ~ 3 • IO^^Mq/H. At 95% of confidence level the 
remnants are hosted by a halo of mass less than 3.6 • IO^^Mq/H. For comparison, the most 
massive halo in the simulation has a mass of 4.4 • IO^^Mq/Zi and there are about 15000 halos 
more massive than 3 • 10^^Mq//i. This is a consequence of the poor correlation between first 
PopIII halos and most massive halos at low redshift. 

Finally, combining the results from all our three simulation boxes, we construct in fig. [6] 
the PopIII star formation rate at high z. The total number of PopIII halos that virialize is 
increasing with redshift, reaching a number density of ~ 0.1{Mpc/h)~^ at ^ ~ 30. In our 
small box simulation, this means that the average density of PopIII halos is ~ 10^'^ per grid 
cell. Therefore there is a very small probability of having two collapsed halos within the 
same cell, an event that would not be captured in our model. 



4. When does the first stars epoch end? 

In Sec. [3] we show that the most massive halos at z = 6 have first light progenitors 
that have been formed when already several thousands of other PopIII stars existed. Are 
these progenitors still entitled to be called first starsl That is, when does the first stars 
epoch end? Here we review the question adopting two different definitions to characterize 
the transition from the first to the second generation of stars, namely (i) a threshold for 
the transition given by the destruction of molecular hydrogen and (ii) a metallicity based 
threshold. 



4.1. Molecular Hydrogen destruction 



One criterion for the end of the first light epoch can be based on the destruction of 
Molecular Hydrogen in the ISM due to photons in the Lyman- Werner ([11.15 : 13.6]e\^) 
energy range emitted by PopIII st ars, if 2 is in fact neede d for cooling of the gas in dark 



matter halos of mass ~ 10^ (e.g. see lBromm fc Larsorul2004j ). The fiux in the Lyman- Werner 
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band is about 7.5% of the ionizing flux (i.e. with an energy range above 13.6 eV). A PopIII 
star is expected to emit a total of about 7.6 • 10^^ photons per solar mass (IStiavelli et al. 
20041 ). so if we assume 3OOM0 as a typical mass we have about 1.7 ■ 10^^ if2-destroying 
photons emitted over the stellar lifetime. Only a fraction ^ 0.15 of these photons can 
effectively destroy an H2, molecule, as the most probable outcome of absorption of a Lyman 
Werner photon is a first decay to a highly excited a vibrational level that later returns 
to the fundamental state, with resulting re-ern itted photons below the 11.15 eV threshold 
JshuU fc BeckwithlllQsi Iciover fc BrandlboOlh . Therefore we estimate that ^ 2.5 • 10^^ jj^ 
molecules will be destroyed by a PopIII star. Given the neutral hydrogen number density 
6. 2-10^^ Mpc~^ this means that a PopIII star destroys H2 over a volume 4-10~^Mj9c^/^, where 
^ is the ratio of molecular t o atomic hydr ogen. Assuming a primordial molecular hydrogen 
fraction ^ ^ 10~^ (e.g. see IPeebleslEoOsh . we obtain that a PopIII star has the energy to 
destroy primordial H2 in a volume of ~ lAj Mpc/h)^ . This number is in bro ad agreement 



with detailed radiative transfer simulations by lJohnson. Greif fc Bromnj (120071 ). From Fig.[6l 
it is immediate to see that by 2; ~ 30 the PopIII number density has reached the critical 
level of 0.1{Mpc/h)~^ and therefore around that epoch the radiation background destroys 
all the primordial Once all the primordial H2 has been cleared the universe becomes 
transparent in the Lyman Werner bands and the new H2 formed during the collapse of gas 
clouds is dissociated by the background radiation. In fact, assuming that the abundance of 



H2 formed during collapse is ^coii ~ 5 • 10"^ (e.g. see iHaiman et al.ll2000l ). this means that 
a collapsing IO^Mq halo produces about 1.3 • 10^^ H2 molecules, a negligible number with 
respect to the 2.5 • 10^^ H2 that are destroyed. Our simple estimate therefore suggests that 
around 2; ~ 30 the star formation rate of PopIII stars in IO^Mq halos is greatly suppressed 
and proceeds in a self-regulated fashion where only a fraction ^ 10~^ of the collapsing halos 
are actually able to cool and lead to the formation of massive PopIII stars. Eventually 
the Lyman Werner background is maintained by PopIII stars formed in more massive halos 
(M ^ IO^Mq), cooled by atomic hydrogen, and, at later times, by PopII stars. 

Inspired by these ideas we set the end of the primordial epoch for PopIII formation at 
the point where the primordial H2 has been destroyed, that is around z ~ 30. Of course 
this is only an order of magnitude estimate and to fully address the feedback due to photo- 
dissociating Lyman Werner photons realistic radiative tra nsfer cosmological simulations are 
needed, which may led even to positive feedback (e.g. see lRicotti et al.ll200ll ). In particular 
our estimate does not take into account the effects of self-shielding and the fact that the 
formation timescale for H2 during the halo collapse may be faster than the timescale for 
photo-dissociation by the background radiation. Thus it is possible that the PopIII star 
formation rate at 2; < 30 is not suppressed as much as predicted by our argu ment. However 



our estimate seems to be in broad agreement with the more realistic model by iHaiman et al. 
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( |2000| ) that predicts the onset of a significant negative feedback at z > 25, depending on the 
assumed efficiency of Lyman Werner photon production. 



4.2. Metal enrichment 



Another possibihty to end the first star epoch is based on a I SM metalhcity threshold . 



However, in this case a clear transition epoch is missing (e.g. see iScannapieco et al.l 12003 
Furlanetto fc Loebll2005l ). This is because metal enrichment, driven by stellar winds whose 
typical velocities are many orders of magnitude lower than the speed of light, is mainly a 
local process. Therefore pockets of primordial gas may exists in regions of space that have 
experienced relatively little star formation, such as voids, even when the average metallicity 
in the universe is above the critical threshold assumed to define the end of the PopIII era. 

In any case this definition provides a longer duration for the first star era. In fact to 
enrich the local metallicity above th e Z = lO'^Zfn threshol d, relevant for stopping PopIII 



formation by chemical feedback (see lBromm fc Larsonll2004l ). one SOOM© SN must explode 
for every ?a 2 • IO^Mq total mass volume (DM+barions), assuming on average a PopIII mass 
of 3OOM0 with yield 0.2. For a Milky Way like halo, this means that about 3000 first stars 
SN are needed to enrich the IGM to the critical metallicity. Accordingly to this definition, 
the Pop III epoch would end within a significant fraction of the total simulated volume 
around z ~ 20, when there is th e collapse of dark ma tter halos originated from 3a peaks at 



the IO^Mq mass scale (e.g., see iMadau fc Reesll200ll ). if the suppression in the PopIII star 



formation rate due to lack of H2 cooling is neglected. A further caveat is that very massive 
stars may end up di rectly in Intermediate Mass Black Holes wi thout releasing the produced 
metals in the IGM jHeeer fc WooslevI I2OO2I : ISantos et al1l2002h . 



5. Growth of the PopIII black hole seeds 

From our investigation it is clear that, before the first PopIII progenitor of the most 
massive halo at z = 6 is born, several thousands of intermediate mass {mBH ~ lO^M©) 
black hole seeds are planted by PopIII stars formed in a cosmic volume that will on average 
host a bright z = Q quasar. This result does not allow to establish an immediate correlation 
between the very first PopIII stars created in the universe and the bright z = 6 quasars, but 
neither does it exclude such a link, as the formation epoch for the quasar seed is still at very 
high redshift {z > 40), when radiative feedback from other PopIII stars already formed is 
unlikely to affect the formation and evolution of the seed (see Sec. 14.11) . Here we investigate 
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with a simple merger tree code what is the fate of the black holes seeds formed up to the 
formation time of the quasar seed and what are the implications for the observed quasar 
luminosity function. 

We assume Eddington accretion for the BH seeds, so that the evolution of the BH mass 
is given by: 

rriBH = mo exp [{t - to) A.az], (4) 



where mo is the mass at formation time to and t^ai is the Salpeter time (ISalpeterlll964j ): 



t^^^ = / "^^^ = 4.507 ■ lOV^^, (5) 
il-e)LEdd 

where e is the radiative efficiency. 

Using Eq. [S] we can immediately see that a difference of At < 2 • lO^yr, that is of Az 10 
at z = 40, in the formation epoch of the BH seed of a bright quasar is not too important in 
terms of the final mass that can be accreted by z = 6, as this corresponds to about half a 
folding tirae. As suming e = 0.1 until z = 6.4, the highest redshift in the SDSS quasar sample 



(jPan et al.l 120041 ). we obtain a ratio of final to initial mass ms/^/mo = 2.62 • 10^ for z = 50 
and mBn/mQ = 1.78 • 10^ for z = 40. Therefore in both cases there has been enough time 
to build up a 2; ~ 6 supermassive black hole with mass rriBH ^ 10^ starting from a PopIII 
remnant. 

This estimate however highlights that only a minor fraction of the PopIII BH seeds 
formed before z = 40 can accrue mass with high efficiency, otherwise the number density of 
supermassive black holes at low redshift would greatly exceed the observational constraints. 
The first BH seeds in the box are distant from each other, so they evolve in relative iso- 
lation, without possibly merging among themselves. Therefore other mechanisms must be 
responsible for quenching accretion of the first BH seeds. Interestingly if we were to assume 
that accretion periods are Poisson distributed in time for each seed, we would not be able 
to explain the observed power law distribution of BH masses at 2; < 6 around the high mass 
end. A Poisson distribution would in fact give too little scatter around the median value 
and a sharp (faster than exponential) decay of the displacements from the mean accreted 
mass. An exponential distribution of the accretion efficiency is instead required to match the 
observed BH mass function. In addition, it is necessary to assume that the duty cycle of the 
BH accretion is roughly proportional to the mass of the halo it resides in. This is sensible, 
since an accretion model unrelated to the hosting halo mass may lead to the unphysical 
result of possibly accruing more mass than the total baryon mass available in that halo. In 
fact, a BH seed formed at z ^ AO is within a halo of median mass ~ IO^^Mq/H at z = 6 and 
a few percent of the seeds may be in halos with mass below ~ IO^Mq/H at that redshift. 
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To explore this possibility we follow t he merging; his t ory of PopIII halos formed at z = 40 
by means of a merger-tree code based on lLacey &: Cold (Il993l ). We implement a BH growth 
based on Eq. \5[ but at each step of the tree we limit the BH mass to itibh ^ V ^tar, where 
rriiiar is the total baryon mass of the halo that hosts the BH. The results are reported in fig. [71 
If the BH growth is not constrained (or only mildly constrained), then a significant fraction 
of the seeds grows above 10^^ Mq, which would result in an unreal istic number density of 
super massive black holes at z = 6. However, if rj r^i 6 ■ 10~^ (like in lYoo fc Miralda-Escude 
2004 : see also lWyithe fc Loebll2003l ). then we obtain an expected mass for the BH powering 
bright z = Q quasar o f ~ 5 • lO^Mf;^, w' hich is in agreement with the observational constraints 
from SDSS quasars (IFan et al.l 120041 ). By fitting a power law function to the BH mass 



function in the range [0.055 : 0.2] • 10^*^ M0 we obtain a slope a ~ —2.6, while the slope is 
a ~ —3.7 in the mass range [0.2 : 1.0] • 10^°Mq, a value that is consistent within the la 
error bar \ yith t he slope of the bright end of the quasar luminosity function measured by 
Fan et al.l tooi ). 



Another contribution to ease an overproduction of bright quasars may be given by the 
suppression of the early growth of the PopIII BH seeds for the first ~ lO^yr after formation, 
that is for about 2tsai ( 1 Johnson fc Bromnj 120071 ). In fact the radiation from a PopIII star 
may evacuate most of the gas from its host halo, so that the subsequent BH growth is 
quen ched until a merger provi des a new gas reservoir to enable growth at near Eddington 
rate ( j Johnson fc Bromml 120071 ). Also the BH seeds situated in more massive halos would 
probably be more likely to replentish their gas supply earlier. 



6. Conclusion 

In this paper we investigate the link between the first PopIII halos collapsed in a sim- 
ulation box and the most massive structures at z ~ 6, with the aim of establishing the 
relationship between the first intermediate mass black holes created in the universe and 
the super-massive black holes that power the emission of bright z = 6 quasars. We show 
that almost no correlation is present between the sites of formation of the first few hundred 
10^ Mq/H halos and the most massive halos at 2; < 6 when the simulation box has an edge of 
several hundred Mpc. Here the PopIII progenitors (halos of mass Mfs ~ IO^Mq) of massive 
halos at z < 6 formed from density peaks that are ~ 1.5a (Mfs) more common than that 
of the first PopIII star in the {512Mpc/h)^ simulation box. These halos virialize around 
Zni ~ 40, to be compared with Zni ^ 48 of the first PopIII halo. 

This result has important consequences. We show that, if bright quasars and super- 
massive black holes live in the most massive halos aX z ^ 6, then their progenitors at the 
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IO^Mq mass scale are well within the PopIII era, regardless of the PopIII termination mech- 
anism. On the other hand, if the itlbh/o' relationship is already in place at 2; = 6, then 
bright quasars are not linked to the remnants of the very first intermediate mass black holes 
(IMBHs) born in the universe, as their IMBH progenitors form when already several thou- 
sands of PopIII stars have been created within the typical volume that hosts a bright z = 6 
quasar. The IMBH seeds planted by this very first PopIII stars have sufficient time to grow 
up to rriBH € [0-2 : 1] • lO^^M© by 2; = 6 if we assume Eddington accretion with radiative 
efficiency e < 0.1. Instead, quenching of the BH accretion is required for the seeds of those 
PopIII stars that will not end up in massive halos at 2; = 6, otherwise the number density 
of supermassive black holes would greatly exceed the observational constraints. One way to 
obtain growth consistent with the observations is to limit the accreted mass at a fraction 
1] ^ 6 ■ 10~^ of the total baryon halo mass. This gives a slope of the BH mass function 
a = —3.7 in the BH mass range rriBH ^ [0-2 : 1] ■ lO^^M©, which is within the la uncertainty 
of the slope of the b right end of the z = 6 quasar luminosity function {a ~ —3.5) measured 



bv lFan et all (120041 ). 



Another important point highlighted by this study is that rich clusters do not prefer- 
entially host the remnants of the first PopIII stars. In fact the remnants of the first 100 
Pop-Ill stars in our medium sized simulation box (volume of {512Mpc/h)^) end up at 2; = 
on halos that have a median mass of 3 • 10^^ Mq/H. This suggests caution in interpreting the 
results from studies that select a specific volume of the simulation box, like a rich cluster, 
and then progressively refine smaller and smaller r egions with the a im of hunting for the 



first lights formed in the whole simulation (see e.g., iReed et al.ll2005l ). Only by considering 
refinements over the complete volume of the box the rarity and the formation ranking of 
these progenitors can be correctly evaluated. 
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A. Tests for the First Light Monte Carlo Method 

To verify the validity of our Monte Carlo method we perform two main tests. First 
we compare the maximum overdensity at the first light halo mass scale identified over the 
whole simulation box using different grid resolutions, including the analytical expectation 
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(that is equivalent to assume that the whole box is just one cell). The results are reported in 
Fig. [8] and confirm indeed that the method is independent of the grid size, with an excellent 
match between all the probability distributions computed. The figure has been obtained 
by first generating a Gaussian random field with a = 0.1 on a = 64^ grid and adopting 
(T/s = 0.008 and k = 4. Then we progressively bin grid cell values to obtain lower resolution 
versions of the original field. The variance in the low resolution grids scales as 
and the values for ajs and k are correspondingly increased. As a second test, shown in Fig.[9l 
we have realized a constrained low resolution (A^ = 128'^) version of the initial conditions for 
our 5*1 simulation and we have then carried out the run down to 2; = 6. From the snapshot 
at z = 6 we identify the most massive halos in this simulations, verifying that there is a 
good spatial and mass match with the original 512^ run. The redshift distribution of the 
first PopIII progenitor for the most massive halos has been computed using our method and 
compared with that of the original run. The agreement is very good, especially considering 
that the dynamics of the dark matter halos has been followed at a resolution 64 times lower. 
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Fig. 1. — Probability distribution functions p{m) for the maximum m oik random Gaussian 
fluctuations representative of mass scale for halos hosting z = Q QSO candidates (red curve) 
and PopIII stars (black curve). To compute p{m) for QSO hosting halos we have assumed 
a box of 720Mpc/h and a mass scale of Mqh = 4.3 IO^^Mq/Zi which implies k = 180^. 
This curve represents the probability distribution for a sigma peak that at z ~ 6 leads to 
one of the most massive halos in the simulation volume. The curve associated to first light 
perturbations (solid black, with Mfs = IQ^Mq/K) is derived using ats = \/o'm, — o'li 

W J s qh 

and k = 160'^: it represents the probability distribution for the maximum of the sub-grid 
scale fluctuations at the Mfs scale within one cell of the 180'^ volume. The dashed green 
line represents the probability distribution of the density fluctuation associated to the first 
PopIII progenitor of a 5.7a{Mgh) peak, m is given in units of CMg^- 
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Fig. 2. — Left panel: Distribution for the ranking of the collapse epoch for the oldest PopIII 
halo progenitor (with Mfs = 10®Mq//i) of the most massive halo (black line) and averaged 
over the 10 most massive halos (red line) at 2; = 6 in the {720Mpc/h)^ box simulation. The 
cardinality is measured over 600 Monte Carlo realizations. Right panel: Like in the left 
panel but distribution of the collapse time Zni for the oldest PopIII progenitor. The blue line 
represents the collapse redshift of the first PopIII star perturbation, while the dotted green 
line refers to the collapse redshift of the 100th PopIII in the box. 
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Fig. 3. — Probability distribution of the collapse time for the oldest PopIII progenitor of the 
most massive halo at 2; = 6 (upper panel) and for the first PopIII halo formed in the box 
(bottom panel) in the {720Mpc/h)^ box simulation. The curves refer to different masses for 
the PopIII halo: Mj, = 3.4 ■ lO^M^/h (blue, short dashed), M/, = 1.0 • IO^Mq/H (black, 
solid) and M/, = 3.0 ■ 10^ Mq/H (red, long dashed). 
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Fig. 4. — Like in Fig. [2l but for a simulation of a {512Mpc/h)^ volume and using 200 Monte 
Carlo realizations. 
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Fig. 5. — Like in Fig.[2l but for the S2 simulation, with volume of {60Mpc/h)^ and ag = 0.75 
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Fig. 6. — Cumulative comoving number density of virialized PopIII lialos {Mfg = IO^Mq/Ii) 
versus redshift for 2; > 29 as measured by means of our MC code. 




Fig. 7. — Final BH mass at z = 6 for a lOOM^/h seed starting Eddington accretion with 
radiative efficiency e = 0.1 at 2 = 40. The growth of the BH is hmited to a fraction rj of the 
total baryon mass of its halo, as obtained with a merger-tree Monte Carlo code. 



- 28 - 



20 



N =64 



15 



b 

^ 10 



5 



N =32 



N =1 







0.6 



0.65 0.7 0.75 0.8 



0.85 



Fig. 8. — Probability distribution for the maximum density fluctuation at first light scale 
(cTfi) in a synthetic test of the refinement method obtained generating first a gaussian random 
field with a = 0.1 on a A^"^ = 64^ grid and then applying our MC method with a refinement 
factor 4, subgrid fluctuations afs = 0.08 and 400 Monte Carlo realizations. We progressively 
downgrade the resolution of the simulation grid to A^^ = 32 and A^^^ = 16, increasing afs 
following our prescription. The MC sampling results are compared with the theoretical 
probability distribution for afi (dashed line). 




Fig. 9. — Probability distribution for the earliest PopIII star in the two most massive halos 
of the SI run obtained from the full resolution (red dashed line) and from a lower resolution 
constrained realization of the same initial conditions (solid black). 
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Table 1: Simulations Summary 



ID 


Lbox h/Mpc 




1,1/3 


Mfs h/Mo 


Nmc 


SI 


60 


0.9 


5 


8.6 10^ 


100 


S2 


60 


0.75 


5 


8.6 10^ 


100 


Ml 


512 


0.9 


40 


1.0 10^ 


200 


LI 


720 


0.9 


57 


1.0 10^ 


600 



Note. — Summary tabic with the details of the numerical simulations carried out in this paper. The first 
column gives the simulation ID ; Li,ox (second column) is the box size, while the third column contains the 
value of cTg used to generate the initial conditions, kref (fourth column) is the refinement factor from the 
mass of a single particle in the N-body run to the mass assumed for a PopIII star halo (M/^; fifth column). 
The last column (Nmc) reports the number of different Monte Carlo realizations of the PopIII halo density 
field for each N-body realization. 



